home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / dbesy1.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  3.5 KB  |  81 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((nty1 0)
  12.       (xmin 0.0)
  13.       (xsml 0.0)
  14.       (by1cs (make-array 20 :element-type 'double-float))
  15.       (twodpi 0.6366197723675814)
  16.       (first nil))
  17.   (declare (type f2cl-lib:logical first)
  18.            (type (simple-array double-float (20)) by1cs)
  19.            (type double-float twodpi xsml xmin)
  20.            (type f2cl-lib:integer4 nty1))
  21.   (f2cl-lib:fset (f2cl-lib:fref by1cs (1) ((1 20))) 0.03208047100611909)
  22.   (f2cl-lib:fset (f2cl-lib:fref by1cs (2) ((1 20))) 1.2627078974335004)
  23.   (f2cl-lib:fset (f2cl-lib:fref by1cs (3) ((1 20))) 0.0064999618999231745)
  24.   (f2cl-lib:fset (f2cl-lib:fref by1cs (4) ((1 20))) -0.08936164528860505)
  25.   (f2cl-lib:fset (f2cl-lib:fref by1cs (5) ((1 20))) 0.013250881221757098)
  26.   (f2cl-lib:fset (f2cl-lib:fref by1cs (6) ((1 20))) -8.979059119648353e-4)
  27.   (f2cl-lib:fset (f2cl-lib:fref by1cs (7) ((1 20))) 3.647361487958307e-5)
  28.   (f2cl-lib:fset (f2cl-lib:fref by1cs (8) ((1 20))) -1.0013743816660006e-6)
  29.   (f2cl-lib:fset (f2cl-lib:fref by1cs (9) ((1 20))) 1.9945396573901739e-8)
  30.   (f2cl-lib:fset (f2cl-lib:fref by1cs (10) ((1 20))) -3.023065601803382e-10)
  31.   (f2cl-lib:fset (f2cl-lib:fref by1cs (11) ((1 20))) 3.6098781569478117e-12)
  32.   (f2cl-lib:fset (f2cl-lib:fref by1cs (12) ((1 20))) -3.4874882972875826e-14)
  33.   (f2cl-lib:fset (f2cl-lib:fref by1cs (13) ((1 20))) 2.7838789715591766e-16)
  34.   (f2cl-lib:fset (f2cl-lib:fref by1cs (14) ((1 20))) -1.867870968619488e-18)
  35.   (f2cl-lib:fset (f2cl-lib:fref by1cs (15) ((1 20))) 1.0685315339116826e-20)
  36.   (f2cl-lib:fset (f2cl-lib:fref by1cs (16) ((1 20))) -5.274721956684483e-23)
  37.   (f2cl-lib:fset (f2cl-lib:fref by1cs (17) ((1 20))) 2.2701994031556638e-25)
  38.   (f2cl-lib:fset (f2cl-lib:fref by1cs (18) ((1 20))) -8.595390353945233e-28)
  39.   (f2cl-lib:fset (f2cl-lib:fref by1cs (19) ((1 20))) 2.8854043798337947e-30)
  40.   (f2cl-lib:fset (f2cl-lib:fref by1cs (20) ((1 20))) -8.647541138937175e-33)
  41.   (setq first f2cl-lib:%true%)
  42.   (defun dbesy1 (x)
  43.     (declare (type double-float x))
  44.     (prog ((ampl 0.0) (theta 0.0) (y 0.0) (dbesy1 0.0))
  45.       (declare (type double-float dbesy1 y theta ampl))
  46.       (cond
  47.        (first
  48.         (setf nty1
  49.                 (initds by1cs 20
  50.                  (* 0.1f0 (f2cl-lib:freal (f2cl-lib:d1mach 3)))))
  51.         (setf xmin
  52.                 (* 1.571
  53.                    (exp
  54.                     (+
  55.                      (max (f2cl-lib:flog (f2cl-lib:d1mach 1))
  56.                           (- (f2cl-lib:flog (f2cl-lib:d1mach 2))))
  57.                      0.01))))
  58.         (setf xsml (f2cl-lib:fsqrt (* 4.0 (f2cl-lib:d1mach 3))))))
  59.       (setf first f2cl-lib:%false%)
  60.       (if (<= x 0.0) (xermsg "SLATEC" "DBESY1" "X IS ZERO OR NEGATIVE" 1 2))
  61.       (if (> x 4.0) (go label20))
  62.       (if (< x xmin) (xermsg "SLATEC" "DBESY1" "X SO SMALL Y1 OVERFLOWS" 3 2))
  63.       (setf y 0.0)
  64.       (if (> x xsml) (setf y (* x x)))
  65.       (setf dbesy1
  66.               (+ (* twodpi (f2cl-lib:flog (* 0.5 x)) (dbesj1 x))
  67.                  (/ (+ 0.5 (dcsevl (- (* 0.125 y) 1.0) by1cs nty1)) x)))
  68.       (go end_label)
  69.      label20
  70.       (multiple-value-bind
  71.           (var-0 var-1 var-2)
  72.           (d9b1mp x ampl theta)
  73.         (declare (ignore var-0))
  74.         (setf ampl var-1)
  75.         (setf theta var-2))
  76.       (setf dbesy1 (* ampl (sin theta)))
  77.       (go end_label)
  78.      end_label
  79.       (return (values dbesy1 nil)))))
  80.  
  81.